/*
 * The main file for computing (timing?) HNF and SNF.
 */

#include "utils.h"
#include <time.h>
#include <fstream>
#include <iostream>

using namespace std;

/*
 * -- Sequential --
 * Algorithm 2 for computing the TBM SNF 
 * http://www.informatik.uni-kiel.de/~gej/publ/parsmith2.pdf
 */
int tbm_seq(int *A, unsigned int n)
{
  int r; // rank
  // Put in HNF

  // IF r<m
   if (r < n)
     {
       /*  // THEN FOR l = 1,··· ,r
       for (int l = 0; l <= r; r++)
	 {
	  // IFl≠ il
	   if (l != A[l]) // Is that right?
	     // THEN al ↔ ail (swap rows)
	     for(int i = 0; i < n; i++)
	       // swapping the row contents l and isubl??
	       ;
	       }*/
       // FOR l = r,··· ,1 (backwards)
       for (int l = r; r > 0; r++)
	 {
	   // FOR s = r + 1, · · · , m
	   for (int s = r + 1; s < n + 1; s++)
	     // COL-ONE-GCD (A, l, l, s)
	     utils::gcdCol(A, l, l, s, n);
	 }
     }
   // FORl=1,···,r−1
   for (int l = 1; l < r; l++)
     {
       // FOR s = r − 1, · · · , l + 1 (backwards)
       for (int s = r-1; s > l; s--)
	 {
	   //  COL-ONE-GCD (A, s, l, s + 1
	   utils::gcdCol(A, s, l, s+1, n);
	   //  ROW-ONE-GCD (A, s, s, s + 1)
	   utils::gcdRow(A, s, s, s+1, n);
	 }
       //  COL-ONE-GCD (A, l, l, l + 1)
       utils::gcdCol(A, l, l, l+1, n);
     }
   utils::bandToDiag(A, n);
   utils::diagToSmith(A, n);
   return 0;
}

/*
 * -- Parallel --
 * Algorithm 2 for computing the TBM SNF 
 * http://www.informatik.uni-kiel.de/~gej/publ/parsmith2.pdf
 */
int tbm_par(int *A, unsigned int n)
{
  return 0;
}

/*
 * Main method for SNF
 */
int main ( int argc, char *argv[] )
{
  if ( argc != 3 ) // argc should be 2 for correct execution
    // We print argv[0] assuming it is the program name
    cout<<"usage: "<< argv[0] <<" <filename> <dimension>\n";
  else {

    // argv[1] is a filename to open
    FILE* the_file = fopen(argv[1], "r+");

    // check to see if file opening succeeded
    if (the_file == NULL)
      perror("Error opening file");
    else
      {
      // initialize dimension and matrix
	unsigned int n = (unsigned)atoi(argv[2]);
	int* A;

        for(int i = 0; i < n; i++)
	  for(int j = 0; j < n; j++)
	    fscanf(the_file, "%i", &A[RIDX(i,j,n)]);

	// explicitly close the file
	fclose(the_file);

	// print out matrix as a test
	utils::printMatrix(A, n);

	hrtime_t seqstart, seqend, parstart, parend, seqtime, partime;
	float savg, pavg;

	// Time sequential execution
	seqstart = gethrtime();
	savg = tbm_seq(A, n); 	  // Triangular Band Matrix Algorithm - sequential
	seqend = gethrtime();
	seqtime = seqend - seqstart;

	// Time parallel execution
	parstart = gethrtime(); 
	pavg = tbm_par(A, n); // Triangular Band Matrix Algorithm - parallel
	parend = gethrtime();
	partime = parend - parstart;

	// Do results match?
	//checkresults(savg,pavg);

	// Print timing results
	double spdup = ((double)seqtime)/((double)partime);
	printf("Seq Time : %lld, Par Time : %lld, Speedup : %f\n", seqtime, partime,
	       spdup);
	return 0;
      }

  }
}
